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ABSTRACT 



. Solar flares are large-scale releases of energy in the solar atmosphere, which are 

, characterised by rapid changes in the hydrodynamic properties of plasma from the pho- 

' tosphere to the corona. Solar physicists have typically attempted to understand these 

I ' complex events using a combination of theoretical models and observational data. From 

o ■ 
in 

c/3 ■ statistically significant comparisons between theory and observations, due primarily to 

^ the large number of free parameters associated with physical models. This class of 
ill-posed statistical pr oblem is ideal l y suit ed to Bayesian methods. In this paper, the 



m . . . 

^ , solar flare studied by IRaftery et al.l ()2009l ) is reanalysed using a Bayesian framework. 

I This enables us to study the evolution of the flare's temperature, emission measure and 

fSJ I energy loss in a statistically self-consistent manner. The Bayesian-based model selec- 

tion techniques imply that no decision can be made regarding which of the conductive 
■ or non-thermal beam heating play the most important role in heating the flare plasma 

during the impulsive phase of this event. 



Subject headings: Sun: corona — Sun: flares — methods: statistical 



1. Introduction 

Solar flares are triggered by the instability of the magnetic field. This can result in the 
direct or indir ect heating of c hromospheric plasma leading to the process known as "chromospheric 
evaporation" (jMilligan et al. 1 [2006a, tbi). There are two ways to provoke evaporation (chromospheric 
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plasma upflow): usi ng thermal energy a nd /or using non-thermal energy according to the energy 
release mechanisms (jMariska et al.lll993l ). In the thermal energy model, heat is unleashed in the 
coronal portion of the loop (possibly because of magnetic reconnection) . Then the thermal energy 
is carried downward to the upper chromosphere via conduction, where the deposited energy heats 
the plasma and stimulates it to move slowly upward. In the non-thermal energy model, the energy 
release (again in the corona) is in the form of a non-thermal electron beam. Electrons then collide 
with dense chromospheric plasma, which they heat, causing the plasma to expand both upward at 
very high speeds and downward at lower speeds. 



This paper will focus on the statistical analysis of the observations outlined in iRafterv et al 



( 20091 ) and will b e compared with resul ts from the Enth alpy Based Th e rmal Evolution of Loops 
(EBTEL) model (jKlimchuk et al.1 boosl ). Comparing to iRaftery et al.1 (j2009l ) this research will 
handle the statistics with care to derive which of the thermal or non-thermal heat flux is more 



dominant, during a sing l e flare event, using different mode l comparison techniques (IKass 



Kass fc Greenhousel[l989l : IXass fc E,aft^ll99d : lR,aftervl[l99fil : Icivde et allboOTl : Uitkin et al 



1993 



20091 ). 



Furthermore, the method used here treats the uncertainties for the time constraints with equal 
respect as the uncertainties of the temperature and emission measure, which are difficult to incor- 
porate into our analysis regarding Classical statistics. The thermal heat flux is proportional to the 
non-thermal heat flux plus a constant background heat flux. The parameter set is extended with 
the addition of the radius of the loop in order to make the model more realistic. Last but not least, 
the ratios between the radiative loss rate of the transition and the corona, the average coronal and 
apex temperature, and the coronal base and apex temperature are assumed to be free parameters 
so that their values will be decided by the data. For more information about these parameters refer 
to Section [3l 

One advantage of the EBTEL model is that we can incorporate both a direct and a non-thermal 
heat input, which is not restricted by the data. Thus, in order to produce a temperature, density 
and/or emission measure profile, a specific form for the non-thermal heat fiux function must be 
assumed. There appears to be a connection between the heating function form and the temperature 
profile, more so than with the emission measure profile. For example, a sudden increase in the heat 
input will lead to an abrupt uplift in the temperature, whereas a flat heating function will lead to 
a smoother temperature proflle. Hence, it is important to understand and evaluate the impact any 
changes in the thermal and non-therm al heat flux have on the plasma evolution. Several forms of 
thermal heat fiuxes have been tested in lAdamakid (|2009l ). Here we present only two of them: Half 
Gaussian profile and Pull Gaussian profile. The Half Gaussian profile represents a sudden switching 
off of an electron beam, whereas the Full Gaussian profile imita tes a gradual one. T he rest of the 
paper has been structured as follows: Section [2] comments on the lRafterv et al.l (|2009l ) observations 
and some assumptions in the EBTEL model. Section [3] discusses the details of the analysis regarding 
the data distribution and the parameter set. Section H] addresses to the choice of the parameters in 
the prior distributions. Section [5] briefiy presents the model selection methods applied here. Section 
El presents results of different models and parameter estimations. Finally, Section [7] summarise the 



- 3 - 



findings from this study and presents how this work can be further progressed. 



2. Previous Work 



2.1. Observations of a C-Class Solar Flare 



The temporal evolution of temperature and emission measure in a C3.0 solar flare observed 
on March 26, 2002 has been analysed. During a typical solar flare, the temperature rises from 
less than one MK up to ~ 20 M K, and then co o ls dow n to the pre-flare temperatures. In the 
example under consideration here, iRaftery et al.l (120091 ) used different instruments to track this 
evolution: the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; > 5 MK), 
GOES-12 (5 - 30 MK), the Transition Region and Coronal Explorer (TRACE 171 A ; 1 MK), 
and the Coronal Diagnostic Spectrometer (CDS; ~ 0.03 — 8 MK). A notation should be made 
that TRACE data were not included in the emission measure analysis. This is as a result of the 
instrument being sensitive to multip le emission lin e s, mak ing it complex to define the contribution 
function. The reader is referred to iRafterv et al.l (|2009l ) which outlines the reasons why and at 
which particular time of the flare each particular instrument was employed. Figure [T] depicts a 
visual display of the observation results together with the underlying uncertainties. 

The flare began with evidence of pre-flare heating at its onset. This was followed by explosive 
chromospheric evaporation during the impulsive phase and gentle chromospheric evaporation dur- 
ing the early decay phase. It is believed that the plasma reached a peak temperature of just more 
than 13 MK in approximately 10 min. Conduction losses dominated over radiative losses for the 
i nitial ^ 300 s of t he post-flare decay, whereas for the next ~ 4000 s, radiative losses prevailed. 



Raftery et al.l (|2009l ) concluded that approximately equal direct and non-thermal heating mecha- 
nisms produced the data observed according only to an estimation of the ratio between the thermal 
and non-thermal heating rate, which is of course very critical. The best fit and parameter intervals 
were derived using an acceptable fit to the data by eye. 



2.2. The EBTEL Model 



The Enthalpy-Based Thermal Evolution of Loops (EBTEL) model OKlimchuk et al.ll2008l ) is 
a zero-dimensional hydrodynamic (OD HD) model that takes into account the significance of in- 
troducing the enthalpy into the system. The difference with the ID HD model is that the ID HD 
model provides us with temperature, density, and velocity profiles along a given magnetic field line 
as time evolves (that is temperature, density, and velocity are a function of both space and time), 
whereas OD HD models describe the average evolution of these values along a coronal strand as a 
function of time only. Subsequently, OD HD models require less computing time than ID models, 
but at the expense of losing the spatial resolution. In spite of this, the EBTEL code has been used 
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Fig. 1. — : Temper ature and emission measure observations of a C3.0 solar flare reported in 



Rafterv et alj ((20091). Top: The EBTEL te mperature evolution that best reproduced the ob- 
servations, according to iRafterv et al.l (|2009l ). Different instruments were employed to measure 
the temperature at different times. Bottom: The EBTEL model and observed emission measure 
evolution. TRACE observations were not included. 
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in a wide range of studies regarding different heat input forms (see iKlimchuk et al.ll2008l . for more 
details). The importance of OD models compared to the ID models is that they can give similar 
results for the plasma response to a sudden heat input, despite the fact that they use up to four 
orders of magnitude less computing time. 

To reflect the physical processes taking place within the standard flare model, the EBTEL 
model allows for both thermal and non-thermal heating of the plasrna in the system. It has 
been noted that proton precipitation is not included (see lAschwandenl |2005| ) . Proton beams are 
expected to excite strong kinetic Alfven waves. The turbulence caused by kinetic Alfven waves 
contains enough energy to produce the non-thermal velocities observed in flares, and therefore 
could contribute to impulsive primary plasma heating between the reconnection regions and the 
flare footpoints. 



2008 ): 



The governing equations of the EBTEL model for the coronal part of the loop are ([Klimchuk et al 
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where P is the pressure, t the time, Q the direct (thermal) heating rate, n the electron num- 
ber density, T the temperature, A(-) the optically thin radiative loss function, the non-thermal 
energy flux (J-" = SJ', with £ the mean energy of the accelerated (non-thermal) electrons imping- 
ing the chromosphere in keV and J the nonthermal particle flux), k the Boltzmann's constant 
(~ 1.38 X 10~^^ erg K~^), L the length of the loop from the coronal base to the apex, Fq the 
thermal conduction (heat flux) at the beginning of the coronal region, Rtr the radiative cooling 
rate at the transition region, ci = the ratio of the radiative loss function between the transition 
region and the corona, C2 = -fr the ratio between the average temperature at the coronal part of 
the loop and the temperature at the apex of the loop, cz = the ratio between the temperature 
at the base of the corona and the temperature at the apex of the loop and the over bars indicate 
spatial averages along the coronal section of the loop. Furthermore, the transition region is treated 
separately. 

Deviations of the observations from the model could be caused by two distinct reasons: First, 
the EBTEL model treats the effects of the non-thermal electron beam in a very simple way. It is 
assumed that all the beam's energy goes into evaporating plasma upwards into the loop, which is 
reasonable for gentle evaporation. However, for explosive evaporation, some of the beam energy will 
be used to drive chromospheric downflows, making the total observed non-thermal energy larger 
than that predicted by the EBTEL model. Second, the flare loop isalmost certainly constructed 
of many strands that are heated at different times. iRafterv et al.l (j2009l ) assume the flaring loop 
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consists of individual strands, all heated simultaneously. Therefore, they model the loop as one 
"fat strand". However, some strands are likely to be heated at a different time compared to this 
solid structure. 



3. Defining the Model Parameters 
3.1. Data Distribution 



The errors in temperature come from the width of the contribution function for the particular 
emission line (for the CDS points) and from the width of the instrument response function for 
GOES and TRACE. The RHESSI temperature error comes from the uncertainty in the fit to the 
Maxwell Boltzmann distribution. 

As there is no information about the form of the errors, it is assumed that the errors between 
the observed temperatures and emission measures and those predicted from the model are normally 
distributed with mean zero and standard deviation derived from the error bars. In this particular 
case study we assume "Sfi" belief, i.e the distance between the mean of the dist ribution and the 



uppe r or lower error bar is three standard deviations of the normal distribution (see lAdamakis et al 



20ld . for a link between error bars and statistical distributions). The reason for choosing a Normal 
distribution as the likelihood function is two-fold: (i) the error bars are symmetric and intuitively 
a symmetric likelihood is needed, and, (ii) it is the standard likelihood distribution that is applied 
in many parametric models in the absence of additional information. 

This means that the combined likelihood which contains both temperature and emission mea- 
sure information is of the form: 

p(D\P) 



(2^) 



exp 



(ni+n2)/2 [j-ini \ „ \ 

^ (Ti - fi(P))' ^ {EMi - EMi{P]^^ 
i=i ^^ii i=i 



2^1 



I (max (T) < 3 X 10^) , 

where D = (T, EM) with T = (Ti, . . . , r„J the temperature data-points and EM = {EMi, EMn^] 
the emission measure data-points, Tj and EMi are the temperatures and emission measures respec- 
tively proposed by the EBTEL model, P is the parameter set (see end of Section[3]), ni is the number 
of observed temperature values, n2 is the number of observed emission measure values, an are the 
standard deviations of the temperature errors, ai2 are the standard deviations of the emission 
measure errors and X(-) is the indicator function which is given by: 



X (max (T) < 3 X lO'^ 



if max (T) < 3 X 10^ K 
otherwise, 
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where max (T) is the maximum temperature value of the coohng curve. The temperature profile is 
restricted to not exceed 30 MK because it is believed that temperatures above this threshold will 
create unphysical data. Temperature and electron density (T and n respectively) are calculated 
numerically from Equations ([I]) - ([3]). First order finite difference methods are employed to solve 
numerically the partial differential equations, after providing initial values for temperature, density 
and pressure. Emission measure (EM) can be computed from 2n^7rr^L, where r is the dimensionless 
radius of the loop and L is the dimensionless length of the loop. 



3.2. Parameters for Non-Thermal Heat Flux 



As discussed in Section [U one of our aims is to compare which of the Half Gaussian and Full 
Gaussian functions for the non-thermal heat flux fits better to the data observed. In the case of 
the Half Gaussian the non-thermal heat flux function is given by the form: 



A 



0, t> 



t < n 



where t is the time, ^ is the time that maximise the non-thermal heat flux function (which can be 
adopted as a parameter), ai is a value that defines the width of this distribution (which can be 
adopted as a parameter as well) and A is the amplitude of the function (parameter). For the total 
non-thermal heat function Ffot we have: 

poo pfi 

J^^^^ = / F{t)dt = / F{t)dt = P{0<t<n)A, 



where P( 



JO JO 

defines the probability. We have P(0 < t < fi) 



p < 0.5 if and only if ai 



where <I>~^(-) is the cumulative distribution function of the Gaussian distribution with mean 
and standard deviation 1. In the case that fi is big and ai is small {e.g. ai < fJ,/3), then 
P (0 < t < /i) ~ ^. For these examples, we can assume J-'tot = 4' Moreover 
dimensionless total non-thermal heat flux, then Ttot = J~tot ^ ^^Ss cm~^. 

In the case of the Full Gaussian the heat flux function is given by the form: 



if Tl^t is the 



A 



27r(Ji 



exp 



(t-/i) 



21 



2al 



We have P{t > 0) = p (with 0.5 < p < 1) if and only if cji 
P{t > 0) ~ 1 and for these cases we can assume J-'tot = A. 



■ If again, e.g., oi < ^/5, then 



3.3. Parameters with the EBTEL Model 



The thermal heating rate is assumed to be of the form: 

Q{t) = aFi{t) + B. 



(4) 
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Here Q is the direct heating rate; B = B' x 10~^ ergs cm~^ s~^ is a constant background heating 
rate that is occurring (B' is dimensionless and is included as a parameter) ; J^i (t) is the non-thermal 
heating rate with J-iit) = where L = L' x 10^ cm is the loop's length from the top of the 
chromosphere to the apex (parameter); and a is a factor that defines which of the two heating 
functions (thermal or non-thermal) is dominant assuming that the background heating rate is 
negligible comparing the other two heating sources. Under this assumption, if a = 1 we have equal 
amounts of thermal and non-thermal heating, if a > 1 then thermal heat flux is more dominant, 
whereas if a < 1 then non-thermal heat flux is more dominant. In particular, we will be examining 
comparisons of the form: Hq : a ^ 1, Hi : a = 1, H2 : a > 1 and H3 : a < 1. Finally, we include the 
dimensionless radius of the loop r' (assuming the loop to be homogeneous), with r = r' x cm, 
as another parameter into our analysis. The importance of r' is that we calculate 2n^7rr^L in order 
to compare it with the observed emission measure values. 

Initially, when this data-set w as firstly assigned for a nalysis, it was suggested that c, should 
be fixed to the values proposed in iKlimchuk et al.l (|2008l ). However, it was interesting to leave 
them as pa rameters to che c k whe ther: i) we will get data profiles with values close to those 
proposed in iKlimchuk et al.l (j2008l ). ii) if not, whether these changes affect other parameters and 
subsequently model comparisons. Therefore, in Section 16. H these values are assumed to be fixed 
to given numbers, whereas in Section 16.21 they are assumed to be free parameters which can be 
determined by the data. 



3.4. Dealing with time 



The data values have uncertainties (error bars) not only on the y axis (temperature and 
emission measure), but on the x axis as well (time). These time errors were calculated as being 
the width of the spline interpolated light cur ves. There haye been some attempts in t he past to 



deal with problems of error bars in both axes ( Winsorl 19461: BerksonI 195ol: Javne^ 



of errors in variable models can be found in ICasella Berge 3 (|l990l l: ICheng fc Van NessI (l99^ ): 



1991 



Reviews 



Fuller! (1198711. For extragal a ctic astronomy applica t ions o f these methods the reader is referred to 
Akritas fc Bershadvl d 19961 ): kellvl (I2OO7I ) : IPatri^ ^200^ . 



m 



A Bayesian solution to the problem using the reversible jump MCMC algorithm can be found 



Henderson et al J 1 120001 ). The latter application is utilised in the current paper and according to 



this, the true timcj should be included as another parameter. 



All in all, the parameter set would be: P = (Pi,P2), with Pi = (ti,...,t„J, where ti is 



^The true value of a variable {e.g. Ttrue) is based on the observed value of that variable {Tabs) with the addition 
of an error (AT), t.e. Ttme = T^bs + AT. 
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the time for the ith observation, and P2 = {L' ,7^1^^, fj.,ai,a, B' ,r') for Section [6TT] and P2 = 
{L' jj^fgf, lj,,ai, a, B' ,r' , €1,02,03) for Section [Ol 



4. Defining the Priors for the Parameters 



A major difference between Classical statistics and Bayesian statistics is that the former does 
not consider any prior information for the parameters, whereas the latter can incorporate any 
available information about t he parameters before we o bserve the data. Further discussion on 
choosing priors can be found in Kass Sz Wasserman J 19961). while the imp ortance of priors in model 
comparison can be found in kassl (I1993I ) and iKass Greenhouse! (| 19891 ) . 

There are three main ways of choosing a prior: 



1. Subjective: the prior expresses the experimenter's personal probability that a parameter lies 
within a specified range. 

2. Objective and informative: the experimenter may have information or historical data prior 
to the experiment being undertaken. 

3. Non-informative: expresses ignorance about the value of a parameter and is usually dominated 
by the likelihood function. 



For the current analysis, we incorporated subjective prior wherever it was possible and non- 
informative in any other case. It is worth noting that in case studies like this, prior distributions 
can prevent parameters from taking unphysical values or values that do not agree with what we 
observe. In this case, prior belief refers to any guess we have about a particular parameter that 
is related to the current observations under consideration, rather than prior belief from previous 
studies. 

The belief we have before we undertake the analysis in Section [6] is summarised in Table [TJ 
We conclude from the observations that L' should be around 3 and probably between 2 and 4. We 
have the option to give more weight to the value of 3 and less weight as we move away from that 
value. Hence, we assume that L' has a Gamma distribution, L' ~ ^(37.64,0.08), where Q{a,/3) is 
a Gamma distribution with shape parameter a, scale parameter /3 and mean a/3, which will give a 
95% probability between 2 and 4 with mode at 3. For the total non-thermal energy J^^^j we do not 
have any information at all, therefore tt (J^iot) 1' where 7r(-) denotes the prior probability density 
function. 

We assume that the flare observations are from the cooling phase of the temperature profile; 
from the highest temperature lines downwards, the temperature drops continually. Thus, it is 
very likely that the main flare energy release will occur at most up until the first observation. 
For this reason we choose the mean of the non-thermal heat flux (fi) to have an upper limit at 
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the lower error bar of the first observation (which is 884 s after the beginning of the observational 
period). Since we do not have any other belief about giving weights to any specific values we assume 
7r(/i) oc I{fj, < 884), where X(-) is the indicator function (see Section [3TT]) . Also, ai ~ ^(6.78, 17.28) 
will give 95% probability for the one standard deviation of the heating function to be between 20 
and 200 with mode at 100. Since we do not not have any knowledge about the background heating, 
we can assume an improper prior, e.g. tt{B') oc 1. In Section [6.21 we use ^(1.70,2.36) as a prior 
for ci and the uniform prior Z^(0, 1) for 02,03. The prior for ci will give ~ 99.90% probability for 
values below 20 and those for C2 and C3 come naturally as they will impose an upper limit of one, 
without favouring any values below unity. Note that the prio r for C2 can be conser vative, as there 
might be evidence that this number is close to the value 0.87 (jKlimchuk et al.ll2008l ). Nevertheless, 
we decided not to include this information as this is not an output from the current observations 
and because we wanted to test whether our simulations will converge to values close to 0.87. 



For the radius of the loop as observed in Figure [T] from iRafterv et al.l (|2009l ) , we assume that 
the ratio R = jj \s probably 1/6. It is almost certainly no more than 1/2; so we give the ratio 
a 99.73% probability to be inside the interval [0,0.5]. For this reason we have R ~ ^(6.03,0.03). 
Since r' = R x L' and R, L' have known prior distributions, we can simulate the prior distribution 
of r' to be ^(5.23,0.12), which will give 95% probability between [0.20, 1.24] with mode at 0.49. 

Since there is no prior information regarding the a parameter, it would be preferable to assign 
an uninformative prior just as with B' . However, due to the fact that it is the "important" 
parameter of the analysis (if we want to test the hypotheses Hi,i = 0, ... ,3), an improper prior 
would lead to Bartlett's paradox ( Lindley 1957 : Bartlett 1957l fl Therefore, for sensitivity reasons, 
results from an informative prior with 50% probability for a > 1 and 50% probability for a < 1 
will be presented as well. For the Hq hypothesis (a 7^ 1), this probability density function will be 
of the form: 

' 0.5, < a < 1 

0.5exp(l — a), a > 1; 



7r(a) 



for the H2 hypothesis (a > 1), 7r(a) = exp(l — a); and for the H3 hypothesis (a < 1), a ~ ZY(0, 1). 
As far as the true time is concerned we can assume that it is normally distributed about the 
observed time with standard deviation obtained from the time error bars. 

Our last assumption involves the initial values of the temperature, density and emission mea- 
sure. For this, we assume the pre-f lare conditions to ha ve an upper limit of 0.5 MK, 6 x 10'' cm~'^ 



and 4 x 10 cm respectively (iRaftery et al.ll2009l ). Since the non-thermal heating rate is al 



most zero in the beginning, then because of Equation the dominant heat input will be the 
background heating rate. Hence, these upper limits can indirectly serve as an upper limit to the 



^In brief, Bartlett's paradox states that the less informative the prior of the "important" parameters is, the more 
the simpler models will be favoured. On the other hand, if the prior of the "important" parameters is very informative, 
then the more complex model will be favoured. This implies that the model selection method can be very sensitive 
to the prior information, which is somewhat bizarre. 



background heating rate. 



5. Model Comparison Methods 
5.1. Direct prior inclusion 



Model comparison methods using Bayesian statistics are reviewed in great detail in lClvde et al 



m 



(120071). An excelle r it des cripti on of how Bay es factor is applied in several problems can be found 



Kass k. Rafteryl (119951 ) and iRafteryl (jl99d ) . The key value for estimating the Bayes factor is to 



calculate the marginal density function: 



p(D|Mfc 



p(D|P)7r(P)dP, 



where is the hypothesis (model) under consideration, p(D|P) is the likelihood function and 
7r(P) is the prior distribution. Roughly speaking, the model with the largest marginal density will 
provide evidence of its favour among a given class of models. Essentially, the marginal distribution 
is nothing more than the belief that we have about the data after we have integrated (averaged) 
over all the parameters, for the specific model under consideration. Then, the Bayes factor for 
testing two hypotheses (M^ and Mi) is given by: 



BFk, 



p(D|Mfc) 
p{T>\Mi) ■ 



(5) 



From Equation ([5]) it is clear that if BFk^i ^ 1 there is more evidence in favour the model, 
if BFf^ i <^ 1 there is more e vidence in favour the M; model, whereas if BF^ i ~ 1 the data do not 
favour either model (but see lKass &: Raftervlll995l . for more information). Thus, the challenge here 
is to calculate the marginal densities. Since they are very difficult to calc ulate analyt ically, we will 

is employed. 



turn to numerical methods. In this paper a two-stage MCMC sampler (jMira 



2001 



Laplace method and importance sampling methods for marginal likelihood estimation as well as 
the MCMC sampler w i th tra nsformation of the parameter set for better convergence are described 
in detail in lAdamakis (2009). 



5.2. Indirect prior inclusion 



Although Bayes factor has an intuitive interpretation, it has two major drawbacks: [i) it is 
heavily dependent on the choice of the prior distribution for the important parameters, and, (ii) 
accurate and efficient computation of the marginal likelihood can b e difficult. Attempts have been 
made to mitigate the first iss ue by using the intrinsic Bayes factor (jBerger &: Pericchilll996l . Il998l ) 
or the fractional Bayes factor (j O ' HaganI 1 1 995l ) . According to these methods, the reasearcher should 
have large enought data, so that the sample can be split into a training set to estimate the posterior 
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distribution of the parameters, which in turn will be used as the prior distribution in the remaining 
sample. However, the Bayes factor still depends to the choice of the training set. For this reason, 
one has to take into account each possible training set and calculate the new Bayes factor as the 
average of all the Bayes factors. The application of this was not feasible for the present study, where 
we have only 13 data-points: 7 temperature measures and 6 emission measures. Therefore, in this 
paper, results from information criteria and posterior deviances will be presented for comparison. 



5.2.1. Information criteria 



The Akaike Information Criterio n (AIC), propo sed by lAkaikd (|1974| ) , and the Bayesian Infor- 
mation Criterion (BIC), proposed bv ISchwarzj (|l978l ) are also employed in the current study. AIC 
propose to choose the model that minimises: 

AIC = — 2(log maximised likelihood) + 2A, 

whereas the latter chooses the model that minimises: 



BIC = — 2(log maximised likelihood) + Alogn, 

where A is the number of the parameters and n is the number of the data-points. The difference 
between two BICs: 



kl 



exp 



\{BICk 



Bid) 



can be viewed as an approximation to the Bayes factor without the researcher's comformable choice 
of priors. However, one has to come to terms with the fact that, indirectly, a prior similar to the 
Jeffreys prior is incorporated (|Kass &: Wassermanlll995l ). The most important difference between 
AIC and BIC is that AIC was designed to find the model that produces estimates of the density 
which is close on avera ge to the true den sit^{§, whereas BIC was designed to find the most probable 
model given the data (jWassermanl l200d ) . 



5.2.2. Posterior deviance 



Another promising method for comparing different models can be to only use t he posterior 



deviance distribution of each competing model (jAitkin et al. 



2009 



Liu AitkinI 120081 ). According 



to this method, suppose that pfc(P^*'|D) is the ith independent MCMC draw of the posterior 
likelihood function given the data D and model k. Furthermore, assume that the tth. independent 
MCMC draw of the posterior deviance function for model k is Dfc(pW |D) = —2 logpfc(p[*] |D), where 
log (•) is the natural logarithm function. From this, the distribution of the difference between two 



■^Close is measured by the KuUback-Leibler distance. 
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model deviances can be calculated as Di^- i = Dj. — Di. Therefore, the posterior deviance difference 
distribution can be used to calculate the probability P[-Dfc,i > /3|D]- For large numbers of /3, the 
higher this probability the stronger evidence there is in favour of model / against model k. As is 
always the case with model selection, subjectivity comes before objectivity. That is, in order to 
objectively choose a model, we need to subjectively specify a rule that will lead us to this choice. 
We present here three different ways that can help a researcher to use posterior deviances for model 
selection. 



1. The choice of the value of (3 can vary for different applications. For example, for (3 = our 
value of interest is translated to the probability of model k having just a higher deviance than 
model I. In order to make an association with Classical statistics, the difference between the 
deviances of two nested models is approximately a chi-square distribution with df degrees 
of freedom, where df = vi — v^, the difference between the number of parameters estimated. 
Then the difference between the two deviances is compared with Xi-a;df-, the value that leaves 
probability a to the right tail of a chi-square distribution with df degrees of freedom. For 
example, if a = 0.0 5 and df = 1 then xo.95;i = 3.84. In this case we can choose (3 = Xi-a;df- 
Aitkin et al.l (120091 ) use /3 = 4.4 as this gives a posterior probability of 0.9 for model I, 



under the assumption of equal prior probabilities on each model. Their suggestion is that if 
P[Dk^i > 4.4|D] > 0.9 then there is quite strong evidence in favour of model I against model 
k. 



2. One can also try to find the value of /3 that gives: 

P[Dk,i > /3|D] = 0.50, /3 > 0. 



(6) 



If we adopt the table in iKass &: Rafteryl (Il995l ) , then we will end up with Table [2] for the 
range of /3. It is worth noting that these numbers are driven more from intuition, rather than 
a scientific justification. 



3. Another way to quantify our beliefs is to calculate the value of 7 that gives 

P[Dk,i > 0|D] =7, 7 > 0.50. 



(7) 



Comparison of the three methods There appears to be a connection between /3,7, d/ and Vk- 
That is, for fixed values of /3 and df the more parameters in a model, the smaller the value of 7, 
in order to obtain the same information. For example, if we adopt the approximation L'fc(PlD) ~ 
D(P) + Xuf,, where D( P) is the frequentist dev iance and Xi^^ has a chi-square distribution with 
z^fc degrees of freedom (jSpiegelhalter et al.ll2002l ). and also f3 = 2,df = 2 then for z^^ = 1 we get 
7 = 0.86, whereas for i/^ = 10, we get 7 = 0.63. Figure [2] depicts how 7 varies with z^^ for 
df = {1,5,10} and /3 = {2,6}. On the other hand, if we bin the values of 7, then /? will vary 
according to 7, (if and z^^. This result implies that in order to quantify our belief one should bin 
the values of /3 or 7 but not both. 
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In Section [6] results for Bayes factors, AIC, BIC and posterior deviances are presented. Regard- 
ing the posterior deviances in particular, although we show results from the above three methods, 
we favour the second method, i.e. comparing the value of /3 that satisfies Equation Q with the 
values in Table [2l 



6. Results 

The aim of this analysis is to address the following questions: 

1. Which of the Half Gaussian and Full Gaussian functions fits better to the data we have? 

2. Which of the thermal and non-thermal heat fluxes is dominant? 

3. Should the q values remain fixed or should they be free parameters? 

The second question is the most important regarding the physics of the system. The other two 
questions are associated more with the statistical analysis. Nevertheless, they can indirectly affect 
which heat fiux is most dominant because they contain relevant parameters. The estimations of the 
marginal densities of the different models that are shown in Figure [3] are presented in Table [3l For 
all the statistical models we have assumed that the mean energy of the accelerated (non-thermal) 
electrons impinging the chromosphere is = 15 keV. Different values of £ have a minor effect on 
the temporal profiles. 

Our decision about a will be based on AIC, BIC and the posterior deviance. Bayes factors are 
presented as well for both proper and improper priors for a. However, since there is no real prior 
information for a, Bayes factors should be used for guidance and not for decisioning. 

6.1. Ci Fixed 
6.1.1. Model selection 



The Ci parameters are fixed at c^ = 4, C2 = 0.87, c.s = 0.72 due to an acceptable overall 



agreement with ID HD simulations (jRafterv et al. 



200S; 



Klimchuk et al.1 [2003) . From Table H and 



with a non-informative prior for a we can derive that the Bayes factor between the Full Gaussian 
and the Half Gaussian heating function is 42.10 using the Importance sampling estimator. This 
provid es "strong" evidence in favour of the Full Gaussian, according to the table in lKass &: Rafterv 
(|l995l ). Similar decision is reached when comparing the marginal densities for informative a prior. 
Comparing the solid and the dashed lines at the left panel of Figure |4] we can conclude that Full 
Gaussian gives a better fit to the data. Also, P[L'5,io > 0|D] = 0.86, P[D5,io > 4.4|D] = 0.64 and 
^[^5,10 > 6.3|D] = 0.50. According to Table [2] there is "strong" evidence in favour of the Full 
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Table 1:: Prior belief about the parameters, 
mode 95% probability prior distribution 





3 


[2,4] 


^(37.64,0.08) 


Hot 


NA 


NA 


oc 1 


M 


NA 


NA 


oc < 884) 




100 


[20,200] 


g(6.78, 17.28) 


a 


NA 


NA 


oc 1 


B' 


NA 


NA 


oc 1 


r' 


0.49 


[0.20,1.24] 


^(5.23,0.12) 


Cl 


NA 


NA 


^(1.70,2.36) 


C2 


NA 


NA 


U{<d,l) 


C3 


NA 


NA 


ZY(0,1) 



Table 2:: Rule of thumb for quantifying our belief when comparing between two models. The values 
of /3 satisfy Equation ([6]). 

/3 Evidence 

0-2 Not worth more than a bare mention 

2-6 Positive 

6-10 Strong 

> 10 Very strong 
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Fig. 2. — : 7 as a function of /3,df and Uf^. If we assume that the values of (3 that satisfy Equation ([6]) 
can be binned as in Table [2] in order to quantify our beliefs, then for a given value of /3 the value 
of 7 that satisfies Equation ^ varies according to df and v^- 




Fig. 3. — : Tree diagram that depicts all the different hypotheses we have used. The purpose of 
using these hypotheses is to compare them and select the one that best describes the data-set we 
have. 
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Table 3:: Model selection criteria according to marginal densities, AIC and BIC produced under 
several hypotheses using the EBTEL model to compare with temperature and emission measure 
profiles (jRafterv et al.ll2009l ). This table produces results for both Cj fixed and free parameters 
for comparison tests. The logarithmic marginal densities are estimated using: 1: Laplace method 
with posterior covariance matrix, 2: Laplace method with robust posterior covariance matrix, 3: 
Importance sampling estimation with the probabilit y density from the first stage of the two-stage 
sampler as the additional probability density. See lAdamakis et al.l (|20ld ) for more information 
about these estimations. 





Log-marginal densities 


AIC 


BIC 


Ci 


Heat 


a prior 


Models 


1 


2 


3 












Ml : 


a / 1 


-796.43 


-797.78 


-797.66 


1589.83 


1597.74 








M2 : 


Q = 1 


-797.38 


-798.19 


-797.37 


1587.22 


1594.57 




HG 


Informative 


Mg : 


Q > 1 


-797.20 


-798.29 


-797.53 


1589.83 


1597.74 








M4 : 


Q < 1 


-797.34 


-798.54 


-797.99 


1590.17 


1598.08 


fixed 




Non-informative 


M5 : 


Q / 1 


-792.88 


-793.68 


-795.95 


1589.83 


1597.74 






Me : 


a / 1 


-793.82 


-795.28 


-794.84 


1582.32 


1590.23 








M7 : 


Q = 1 


-794.45 


-795.39 


-794.40 


1581.13 


1588.48 




FG 


Informative 


Ms : 


a > 1 


-793.87 


-794.93 


-794.42 


1582.32 


1590.23 








Mg : 


a < 1 


-794.17 


-795.31 


-794.73 


1582.73 


1590.64 






Non-informative 


Mio 


: a / 1 


-789.07 


-790.48 


-792.21 


1582.32 


1590.23 








Mil 


: a / 1 


-793.26 


-795.98 


-798.47 


1573.52 


1583.12 








M12 


: a = 1 


-793.36 


-796.26 


-798.06 


1571.88 


1580.92 




HG 


Informative 


Mi3 


: a> 1 


-794.07 


-795.91 


-797.96 


1573.52 


1583.12 








Mi4 


: a <1 


-793.23 


-796.60 


-797.62 


1573.98 


1583.58 






Non-informative 


Mi5 


: a / 1 


-789.65 


-792.36 


-796.98 


1573.52 


1583.12 


free par. 






M16 


: a / 1 


-791.06 


-793.25 


-794.28 


1568.66 


1578.26 






Informative 


Mi7 


: a = 1 


-790.26 


-792.38 


-794.39 


1567.62 


1576.66 




FG 




M18 


: a> 1 


-789.99 


-792.37 


-794.32 


1568.66 


1578.26 








Mi9 


: a <1 


-791.48 


-793.32 


-794.71 


1569.05 


1578.65 






Non-informative 


M20 


: a / 1 


-786.85 


-788.49 


-792.02 


1568.66 


1578.26 
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Gaussian. Furthermore, both mformation criteria prefer the Full Gaussian model. Finally, the first 
column of Figure [5] shows the best fit of the parameters regarding the posterior distributions. Even 
by eye we can distinguish between the Full Gaussian and the Half Gaussian function. 

From all the above, we can conclude that there is enough evidence to support that the Full 
Gaussian function for the heating profiles is much more adequate than the Half Gaussian, at least 
for the particular data-set under analysis. 

In order to distinguish between the thermal and non-thermal heat fluxes, we have calculated 
P[Dj^s > 0.6|D] = 0.50, P[Dg^7 > 0.6|D] = 0.50 and P[Dq^8 > 1-2|D] = 0.50. Otherwise, 
P[Dj^s > 0|D] = 0.54, P[Dgj > 0|D] = 0.54 and P[£>9,8 > 0|D] = 0.58. According to Table E] 
none of the models is "worth more than a bare mention". On the other hand, both information 
criteria choose the equal amounts of thermal and non-thermal heating model (My), although its 
AIC and BIG values are very close to the ones of Mg and Mg. Furthermore, the marginal densities 
of Mt, Ms, Mg are very similar to each other, something that supports the idea that the difference 
between the models is "not worth more than a bare mention" . Finally, according to the quantiles of 
Table m the 95% credible interval with informative prior for a contains 1, whereas the 95% credible 
interval with non-informative prior for a does not contain 1. In other words, P{a > 1) is close to 
0.50 with informative prior for a, whereas it is close to 1 with non-informative prior for a. 

In summary, it is very difficult to distinguish which of the two heating forms is more dominant. 
Better data are required to address this issue. Since the Full Gaussian is more preferable than the 
Half Gaussian we will choose to adopt the Full Gaussian results for estimating the parameters. 



6.1.2. Parameter estimation 
The results of this analysis with non-informative prior for a can be vi ewed in Table [ 5l We 



estimate the mean of L to be 29.20 Mm, which is close to what was used by lRaftery et al.l (120091 ) 



(i.e. 30 Mm). The mean total non-thermal heat flux is calculated to be around 79.83 x 10^ ergs 
cm"'^, while the mean time for the peak of this function is 853.5 s after the beginning of the 
flare. For the a parameter a 95% credible interval gives [1.30,101.30], whereas the probability 
P{a > 1) = 55.86% with informative prior for a. We expected the posterior probability of a > 1 to 
be close to 1/2, as the Bayes factor could not provide enough evidence in favour of the hypotheses 
a > 1 or a < 1 (Table [3]). An interesting point to note here is that if we use non-informative 
prior for a, we can conclude that since a = 20.90 then the a > 1 hypothesis is preferable. This is 
something that is misused in astrophysics and can lead to false conclusions. More discussion about 
this is presented in Section [71 
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Table 4:: Posterior distribution information for the a parameter. Results are split according to 
Ci (fixed/parameters), heating flux (Half Gaussian/Full Gaussian) and a prior (informative/non- 
informative). 





Quantiles 


Pin- > 1) 




IICMl 




2.F/A 


50% 


97.5% 


fixed 


HG 


Informative (Mi) 


0.16 


0.91 


3.11 


0.44 


Non-informative (M5) 


0.34 


8.10 


48.47 


0.90 


FG 


Informative (Mg) 


0.22 


1.10 


3.70 


0.56 


Non-informative (Afio) 


1.30 


5.74 


101.30 


~ 1 


free par. 


HG 


Informative (Mn) 


0.07 


0.90 


2.93 


0.45 


Non-informative (M15) 


0.42 


11.63 


32.31 


0.91 


FG 


Informative (Mie) 


0.13 


1.13 


4.08 


0.57 


Non-informative (M20) 


1.13 


23.48 


47.09 


0.98 



Full Gaussian with c, parameters 



c, fixed / HG 



c, fixed / FG 

C/ parameters / HG 

Ci parameters / FG 



1540 1550 1560 1570 1580 1590 1600 
Deviance 




1550 1555 
Deviance 



1560 1565 



Fig. 4. — : Left: Cumulative distribution functions for posterior deviance for models with Cj fixed 
and Half Gaussian (solid), q fixed and Full Gaussian (dashed), Cj parameters and Half Gaussian 
(dotted), Ci parameters and Full Gaussian (dot-dashed). All curves are with non-informative prior 
for a. The further to the left, the better the model. Right: Cumulative distribution functions for 
posterior deviance for Cj parameters and Full Gaussian are presented for a closer comparison. The 
dashed line (a = 1) over plots the solid line [a / 1). The dotted line (a > 1) provides a slight 
better fit than the dashed line and the dot-dashed line {a < 1) provides a slight worst fit than the 
dashed line. However, all the hypotheses provide very close lines, an indication that distinction 
between these models will be very difficult. 
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Tabic 5:: Summary of the posterior inference for botli q fixed and free parameters with non- 
informative prior for a. Results steam from a Full Gaussian non-thermal heat flux profile. 





mean 


mode 


s.d. 


2.5% 


50% 


97.5% 




L' 


2.92 


2.86 


0.12 


2.69 


2.91 


3.19 




tot 


79.83 


44.72 


62.96 


3.77 


75.73 


179.88 






853.5 


881.1 


31.6 


763.0 


863.8 


883.1 


Ci fixed 




101.8 


94.2 


26.9 


54.4 


99.6 


154.6 




a 


20.90 


8.46 


28.23 


1.30 


5.74 


101.30 




B' 


0.19 


0.25 


0.11 


0.01 


0.20 


0.36 




r< 


0.47 


0.51 


0.12 


0.25 


0.45 


0.73 




L' 


3.16 


2.80 


0.43 


2.43 


3.13 


4.12 




■pi 


31.99 


8.35 


40.11 


4.84 


13.57 


153.98 






864.6 


880.6 


18.7 


813.2 


870.3 


883.5 






92.6 


82.0 


24.2 


51.8 


90.5 


146.2 




a 


23.51 


38.09 


15.94 


1.13 


23.48 


47.09 


Ci free parameters 


B' 


0.12 


0.21 


0.07 


0.01 


0.17 


0.24 




r' 


0.67 


0.57 


0.17 


0.39 


0.66 


1.05 




Cl 


2.07 


1.39 


0.77 


0.81 


1.99 


3.79 




C2 


0.86 


0.88 


0.07 


0.73 


0.87 


0.99 




C3 


0.75 


0.69 


0.13 


0.48 


0.75 


0.98 
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Fig. 5. — : Temperature evolution {first row), emission measure evolution {second row) and tem- 
perature against emission measure {third row) using the EBTEL model that best reproduced the 
observations (maximised the posterior distributions). Solid lines represent the fit from the Full 
Gaussian function for the thermal and non-thermal heat inputs, whereas dashed lines represent 
the fit from the Half Gaussian function. For all the curves a non- informative prior for a was 
employed. The first column depicts solutions using the fixed, whereas the second column depicts 
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6.2. Cj Parameters 

6.2.1. Model selection 

We follow the same procedure as in Section 16.1.11 in order to derive the models that best 
describe the data. Regarding which of the two non-thermal heat fluxes is more dominant, there 
is no question that the Full Gaussian is more preferable as can be seen from Table [3] (the Bayes 
factor in favour of the Full Gaussian is 142.54 with a non- informative prior for a, according to the 
Importance sampling estimator). This can be characterised as "very strong" evidence in favour of 
the Full Gaussian. The same preference for the Full Gaussian function can be derived by comparing 
the log-marginal densities for a ^ \ with an informative prior for a. Looking at the dotted and 
dot-dashed lines at the left panel of Figure H] we can conclude that Full Gaussian gives a better fit 
to the data. Also, P[Di5,20 > 0|D] = 0.84, P[L>i5,2o > 4.4|D] = 0.66 and P[I)i5,20 > 7.3|D] = 0.50. 
Again, this is "strong" evidence in favour of the Full Gaussian. Furthermore, both information 
criteria favour the Full Gaussian model. Finally, even by eye, the Full Gaussian model produces 
better results in the second column of Figure [5l Therefore, the fact that we included the Cj as free 
parameters, did not alter the outcome compared to the analysis in Section [6.1.11 where the Cj are 
fixed. 

Regarding the thermal and non-thermal heat fluxes, from the right panel of Figure H] it is 
apparent that there is a slight better fit to the data with the a > 1 hypothesis. Nevertheless, this 
improvement in fit does not seem to be substantial. We have also calculated P[Z)i7^i8 > 0.5|D] = 
0.50, P[Dig^n > 0.5|D] = 0.50 and P[Di9^is > 1.0|D] = 0.50. Otherwise, P[Dn^is > 0|D] = 0.52, 
P[i:>i9,i7 > 0|D] = 0.52 and P[-Di9,i8 > 0|D] = 0.55. According to Table [2] none of them is "worth 
more than a bare mention". On the other hand, both information criteria favour M17, just as when 
Cj are fixed. Moreover, the marginal densities of Mii, Mis, Mig are very similar to each other, 
something that supports the idea that the difference between the models is "not worth more than 
a bare mention" . Finally, according to the quantiles of Table [U the 95% credible interval with 
informative prior for a contains 1, whereas the 95% credible interval with non- informative prior 
for a does not contain 1. In other words, P{a > 1) is close to 0.50 with informative prior for a, 
whereas it is close to 1 with non-informative prior for a. Therefore, we reach the same conclusion 
as in Section [6.1.11 although it is clear that the Full Gaussian function in preferable, the data are 
not sufficient in order to distinguish between the different heating mechanisms. 

6.2.2. Parameter estimation 

All the estimations from the posterior distributions of the parameters can be viewed in Table 
[5l using the Full Gaussian function. We estimate the mean of L to be 31.6 Mm, the mean of the 
total non-thermal heat flux is 31.99 x 10^ ergs cm~^, the mean time for the peak of this function 
is 864.6 s after the beginning of the flare, the mean of the standard deviation of the Full Gaussian 
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is 92.6 s, the mean of the a parameter is 23.51, the mean background heating rate is 0.12 x 10 
ergs cm~^ s~^ and the radius of the loop is 6.7 Mm. 

Regarding the Cj ratios, the mean of the ratio between the radiative loss rate of the transition 
and the corona (ci) is 2.07, the mean of the ratio between the average coronal temperature and 
the apex temperature (02) is 0.86 and the mean of the ratio between the coronal base temperature 
and the apex temperature (03) is 0.75. Once again, we should bear in mind that although the 
estimation for a is greater than unity (23.51) with non-informative prior, a more detailed analysis 
shown in Section 16.2.11 suggests that we cannot distinguish which of the thermal or non-thermal 
heat fluxes is more dominant. Furthermore, assuming informative prior for a, a is greater than 
unity with probability 57.13%. 

6.3. Comparison Between the Hypotheses: a Fixed and a Free Parameters. 

Considering the non-thermal heat flux, both hypotheses (cj fixed and q free parameters) 
propose the Full Gaussian statistical model. However, none of them can distinguish which of 
the two heating mechanisms is dominant (if any). Regarding Question 3, from Table [3] we reach 
different conclusions depending on the estimation method for the marginal distribution: Importance 
sampling does not seem to favour any hypothesis, but the other two methods favour the hypothesis 
Ci free parameters. This can be concluded when a comparison between models Miq and M20 (or 
Mq and Miq) is made. On the other hand, the posterior deviances of Figure H] depict preference 
to the models where q are set as free parameters and the form of the non-thermal heat flux is 
Full Gaussian. We have also calculated P[I?io,20 > 19.2|D] = 0.50, P[£'io,20 > 0|D] = 0.995 and 
P[DiQ^20 > 4.4|D] = 0.98 which indicates a "very strong" evidence in favour of model M20 against 
model MiQ. Similar conclusions can be driven with the information criteria, as both of them show 
clear preference to the Cj free parameter models. Last but not least, the difference between the two 
hypotheses in Figure [6] does not seem to be very important for the thermal evolution (first row), at 
least by eye. However, for the data-set under consideration, the difference is more profound for the 
emission measure evolution (second row). This indicates that the Cj parameters affect the emission 
measure values more than they affect the temperature values. The fact that there are better values 
for Cj than 4, 0.87 and 0.72 is even more clear in the third row where temperature is plotted against 
emission measure. 

All the above indicate that the values introduced for Cj when they are fixed are not very good, 
in terms of maximising the likelihood function. The information added when Cj are free parameters 
is greater than the price we have to pay for introducing three additional parameters. 
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Fig. 6. — : Temperature evolution {first row), emission measure evolution {second row) and tem- 
perature against emission measure {third row) using the EBTEL model that best reproduced the 
observations. Solid lines represent the fit when q are floating as free parameters, whereas dot- 
dashed lines assume Cj to be fixed numbers. Each curve depicts the best fit that maximised the 
posterior distribution of the Full Gaussian statistical models. A non-informative prior for a was 
employed. Temperature is measured in K, while emission measure is in cm~^. 
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7. Discussion 

In order to reveal the mysteries of the Sun, we can break down our investigations into three 
major stages: (i) constructing theoretical models, (ii) gathering observations, and (iii) applying a 
statistical analysis in order to compare different statistical models and/or to restrict the parameters 
of the models. All of them are of equal importance and give us confidence in our results. 

Making inferences solely from the mean or mode of TableOwould warrant utmost criticism. For 
model comparison purposes the 95% credible interval is more robust. However, even with the 95% 
credible interval we can not compare with the a = 1 hypothesis. Also, we do not directly include 
the likelihood function in our calculations. Model comparison techniques using Bayesian statistics 
{e.g. Bayes factor) take into account both the spread of the parameter posterior distributions 
and the information of the likelihood function. However, Bayes factor is sensitive to the choice of 
prior distributions and great attention should be paid when applying this method. On the other 
hand, techniques that are based only on posterior likelihood distributions do not depend on the 
prior distributions of the parameters — at least not directly. We strongly encourage researchers to 
employ and compare various statistical techniques when embarking on the subtle topic of model 
selection. 

In this paper, temperature and emission measure profiles produced by the EBTEL model were 
compared with solar flare observations. The data distribution, the parameter set and the priors 
employed in this analysis were described in detail. The form of thermal and non-thermal heat 
input is much bett e r des cribed using Full Gaussian energy profile than Half Gaussian, which is 



what iRaftery et al.l (120091 ) also used. 



Apart from choosing which energy profile function was more appropriate for the data-set we 
analysed, we were also interested in determining which of the thermal or non-thermal heat fluxes 
was more dominant. The data obtained were not able to provide an answer with great confidence. 
More data-points may be required in order to address this question. It has also been suggested 
that it might be possible to have mostly non-thermal heating in the impulsive phase and mostly 
thermal heating thereafter. To test this, Equation (jl]) can be replaced by: 



Q{t) 



aiTi{t) + B, t<ti 
a2Ti(t) + B, t>ti, 



where ai,a2 and ti parameters under consideration. If the above statement were true then one 
would expect that ti should be close to where the temperature peak of Figure [5] is observed, as well 
as «! < 1 and 02 > 1. 



The Cj parameters were given fixed values in the iKlimchuk et al.l (|2008l ) paper so that the OD 
HD model will approximate the ID HD model. If we assume that the range of these parameters 
provides sufficient approximations, then we should include them as free parameters. In any case, 
we disagree with fixing some parameters to certain values, in order to reduce the parameter set, as 
this might affect the results — unless we have high confidence about these fixed numbers. 
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The fact that the Bayes factor is not so decisive in choosing between the hypotheses Cj fixed 
or free parameters is partially an outcome of the conservative prior distributions we have chosen 
for the Cj parameters (see Section . If we had added more information in the prior distribution of 
C2, say ^(0.5, 1) or even better the Beta distribution ;S(38.49, 5.75), this would have been in favour 
of the Cj parameters hypothesis. In comparison, posterior likelihood techniques and information 
criteria clearly show that Cj should be set as free parameters. More improved estimations for these 
parameters can be seen from the mean or mode of Table [5] (for the particular data-set we analyse). 

Finally, an obstacle presented in this analysis was that of constraining the profiles produced by 
the EBTEL model. For example, we might not want to restrict the initial values of temperature, 
density and/or emission measure profiles. This can produce model profiles that are closer to the 
data profiles, but the initial values might take exceptionally high numbers. For instance, we had 
undertaken the same analysis without fixing the initial values of temperature, density and emission 
measure. This resulted in model profiles with initial temperaturecl of ~ 3 MK, initial electron 
density of ~ 500 x 10'' cm~^ and initial emission measures of ~ 3 x 10^'' cm~'^. This was because 
the background heating rate was three orders of magnitude higher than that in Table El Apart from 
unrealistic estimations of the parameters of interest, this could have also led to unreliable Bayes 
factor estimations with false conclusions regarding the three posed questions in the beginning of 
Section \6\ Naturally, the quality of the output of the analysis is dependent of the quality of the 
input. 

Apart from improved statistical techniques, of equal importance is that improvement upon the 
observations should be made. This means that future missions with new instrumentation should 
provide data-sets with a large enough number of observations in order to distinguish between 
different heating mechanisms. The data-set under consideration provided information only upon 
the decay phase of the temporal evolution. However, a better time resolution for the rise phase of 
the temperature will be needed in order to provide a better estimate for the form of thermal/non- 
thermal heat fiux. And of course, a large sample of solar fiares will be required. 

An assumption made in this paper is that the thermal and non-thermal fluxes have the same 
form, based on a lack of information on the thermal distribution. However, it would be interesting 
to test fluxes of different forms that do not depend on each other. Additionally, a further im- 
provement in the EBTEL model is required regarding the non-thermal heat flux, as it is efficient 
for gentle chromospheric evaporation but suffers from inadequately representing explosive chromo- 
spheric evaporation. Last but not least, several other forms of heating input, like proton beams, 
could be included in an attempt to make the model more realistic. 

SA has been supported by a STFC grant. CLR is supported by an ESA/Prodex grant, ad- 
ministered by Enterprise Ireland. The authors would like to express their gratitude to the referee 



^Initial temperature is assumed to be at the base of the transition region. 
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